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ABSTRACT 

Using self-consistent magnetohydrodynamic (MHD) simulations, we explore 
the hypothesis that nonlinear MHD waves dominate the internal dynamics of 
galactic molecular clouds. Our models employ an isothermal equation of state 
and allow for self-gravity. We adopt "slab-symmetry," which permits motions 
vi and fields Bi perpendicular to the mean field, but permits gradients only 
parallel to the mean field. This is the simplest possible geometry that relies on 
waves to inhibit gravitational collapse along the mean field. In our simulations, 
the Alfven speed va exceeds the sound speed c s by a factor 3 — 30, which is 
realistic for molecular clouds. We simulate the free decay of a spectrum of 
Alfven waves, with and without self-gravity. We also perform simulations with 
and without self-gravity that include small-scale stochastic forcing, meant to 
model the mechanical energy input from stellar outflows. 

Our major results are as follows: (1) We confirm that the pressure associated 
with fluctuating transverse fields can inhibit the mean-field collapse of clouds 
that are unstable by Jeans's criterion. Cloud support requires the energy in 
Alfven -like disturbances to remain comparable to the cloud's gravitational 
binding energy. (2) We characterize the turbulent energy spectrum and density 
structure in magnetically-dominated clouds. The perturbed magnetic and 
transverse kinetic energies are nearly in equipartition and far exceed the 
longitudinal kinetic energy. The turbulent spectrum evolves to a power-law 
shape, approximately v k ~ B\ k j Airp oc k~ s with s ~ 2, i.e. approximately 
consistent with a "linewidth-size" relation a v (R) oc R 1 ^ 2 . The simulations 
show large density contrasts, with high density regions confined in part by the 
pressure of the fluctuating magnetic field. (3) We evaluate the input power 
required to offset dissipation through shocks, as a function of c s /va, the velocity 
dispersion a v , and the characteristic scale A of the forcing. In equilibrium, the 
volume dissipation rate is 5.5(c s /fA) 1 ^ 2 (A/L) -1 / 2 x pa^/L, for a cloud of linear 
size L and density p. (4) Somewhat speculatively, we apply our results to a 
"typical" molecular cloud. The mechanical power input required for equilibrium 
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(tens of L ), and the implied star formation efficiency (~ f%), are in rough 
agreement with observations. Because this study is limited to slab symmetry 
and excludes ion-neutral friction, the dissipation rate we calculate probably 
provides a lower limit on the true value. 



1. Introduction 

The internal dynamics of star-forming galactic molecular clouds is not yet understood. 
Two central questions are (1) what prevents the clouds and their subcomponents from 
collapsing under their own weight; and (2) what generates and controls the turbulent 



fluid velocities that broaden molecular lines far beyond the thermal speed c s (e.g. |Slru 



et al. (1987)|) . One model which has been proposed (e. g. |Scalo fc Pumphrey (1982)]) 



is that the clouds are comprised of clumps on essentially ballistic, collisionless orbits. 
However, while clouds are observed to be clumpy, the volume filling factor of clumps in 
the clouds / ~ 0.03 - 0.08 (e.g. |Perault, Falgarone, fc Puget (1985)| ; |Williams, Bllt^ 



Stark (1995)| ) implies a clump-clump collision time t C oiiis < (4/3)-R c i ump /(/f c i ump ) ~ 10 7 yr, 



which makes the clouds at most marginally collisionless over their lifetimes ( |Blitz fc Shu 



(1980)| ). The clumps are not themselves thermally supported, and they appear to have 
larger internal filling factors and smaller ratios of internal collision time to dynamical 
time. Although internal velocities may be generated by a cloud's self-gravity, purely 
hydrodynamic turbulence - either clumpy or smooth - cannot in itself support a structure 
for longer than the effective collision time (equal to the eddy-turnover time for a uniform 



fluid) because it would dissipate in shocks (see [Elmegreen ( 1985")] and references therein). 



The orbiting-clump model therefore probably cannot account for the internal dynamics of 
molecular clouds at all scales. [] Rather than assuming a clumpy mass distribution a priori, 
it seems better to start with a full fluid model with a compressible equation of state, so 
that clumping can be treated self-consistently. Such a model must have some internal stress 
far more powerful than gas pressure in order to control supersonic motions. 

For some time, magnetic fields have been considered the leading candidate for 
mediating clouds' internal motions and counteracting gravity (see the recent reviews of 
|Shu et al. (1987)] ; |McKee et al. (1993)| ). Magnetic processes have also been identified as 



likely instruments for generating density structure within clouds (e.g. [Elmegreen (1990") 



|Elmegreen (1991)1) , which is observed at all scales down to the limiting telescopic resolution 



1 In particular, it is not apparent how a collisionless system could generate non-self-gravitating clumps 
whose internal velocity dispersions scale with their size (cf. |Bcrtoldi fc McKee (1992)] ). 
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( [Falgarone, Phillips, fc Walker ( 1991 )| ; [Falgarone, Puget, fc Perautt ( 1992 )|) . Measured field 



strengths Bn based on OH Zeeman splittings are in the range 10 — 30/iG ( |Crutcher et al 



for the line-of-sight field in moderate-density regions n# 2 ^ 1000cm -3 (for random 
orientations the mean total field strength is twice as large). Fits incorporating additional 
data from weak-field, low-density HI Zeeman splitting and strong-field, high-density OH 
maser Zeeman splitting yield B ~ 1.5(n^ 2 /lcm~ 3 ) 1 / 2 /iG ( Heiles et al. (1993)1 , and references 
therein). Based on these data, the magnetic field has an energy density comparable 
to the kinetic (and gravitational) energy densities, and therefore can be dynamically 
important. More specifically, Myers fc Goodman (1988a)| show that magnetic, kinetic, 



and gravitational energies are comparable in detail for several clouds at a range of scales, 
suggesting virial equilibrium. The field topology within molecular clouds remains uncertain. 
In optical wavelengths, the linear polarization directions of background stars shining 
through low-density regions undulate smoothly across cloud complexes (e.g. |Moneti et al| 



(1984)| ). To trace higher-density gas within clouds, longer wavelengths are needed. Maps 



of polarized lOO/i thermal emission in several high- mass star-forming regions (( Potsori 



1995)|) , |Hildebrand et al (1995)| , Pavidson et al ( 1995 )| ) also show orderly variation across 



the cloud. If in both cases the polarization is caused by field-aligned dust grains, the data 
imply smoothly-varying mean fields. These preliminary indications on field geometry, if 
confirmed, permit a conceptual separation into cloud support perpendicular to, and parallel 
to, a slowly-varying, untangled, mean field. 

To date, most theoretical work on magnetic fields in star-forming regions has 
concentrated on the role of smooth fields in quasi-static equilibria or configurations 
undergoing laminar rotation and/or collapse (see the reviews of |Nakano (1984)| ; 



Mouschovias ( 1991 )| ; McKee et al. (1993)| ). The absence of turbulent velocities tWb 
exceeding c s in the small, dense cloud cores observed to be the sites of low-mass star 
formation (see, e.g. [Fuller fc Myers ( 1992 )| ) makes them amenable to quasistatic theories. To 



the extent that turbulent magnetic and Reynolds stresses can be included via a barotropic 
pressure, such calculations can also be applied to cases where v tm h > c s . Axisymmetric 
calculations of field-frozen equilibria have quantified the importance of field support 
perpendicular to the mean field direction, which can be expressed succinctly in terms 
of the mass-to-magnetic flux ratio, M/<3> QMouschovias (1976a)| ; Mouschovias fc Spitzer 



1976] ; Tomisaka, Ikeuchi, fc Nakamura (1988)| ). The value of this evolutionary invariant 



determines whether or not an equilibrium can be sustained. 

While static or time-averaged fields are likely key to cloud support at both small 
and large scales, they do not oppose gravity in the mean field direction, and by definition 
cannot produce a large velocity dispersion. For clumps within clouds (reviewed by |Blitz 



(1993)1 ; see also Bertoldi fc McKee (1992)1) , and massive cloud cores (e.g. |Caselli fc Myers| 
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(1995)), however, molecular line observations exhibit linewidths in excess of c s . The inferred 
hypersonic bulk velocities were attributed to MHD waves shortly after their discovery 
QArons fc Max (1975)| ). For Alfven waves, the fluctuating component of the field provides a 



pressure that acts along the mean field, and can therefore oppose gravity in that direction 
( Pewar (19701 ; [Shu et al. (19871 Fatuzzo fc Adams (1993)| ; |McKee fc Zweibel ( 1995)0 . 



The theory of Pewar (1970)| calculates the influence of small-amplitude MHD waves on 
the background state of the fluid, using a locally-averaged Lagrangian. For Alfven waves, 
the effect of the waves is expressed through an isotropic wave pressure P wave = { |<5B | 2 ) /87r, 



where the magnetic disturbance is 5B. Recently, |McKee fc Zweibel (1995)] have used 



Dewar's theory to show that small-amplitude Alfven waves propagating along a density 
gradient obey P wa , vc oc p 1 / 2 (implying effective polytropic index 7 p = 1/2), while waves 
trapped in a contracting cloud obey P wav e oc p 3 ^ 2 (implying effective adiabatic index 
7 a = 3/2). Since gas spheres are dynamically stable to adiabatic perturbations when 
7 a > 4/3 (e.g. |Cox ( 19801) ? large amplitude Alfven waves could potentially support a cloud 
against collapse if they suffered minimal decay, or loss (cf. |Elmegreen ( 19851) into the 
surrounding medium, and obeyed the same scaling. A crucial unknown is the decay rate of 
arbitrary amplitude MHD waves in conditions appropriate for a molecular cloud. If Alfven 
waves are responsible for the internal linewidths of molecular clouds and support against 
gravity and external pressure, then any decay must be replenished if a quasi-equilibrium 
is to be maintained. Prior to high-mass star formation, the ultimate source of new wave 
energy must be the gravitational potential of the cloud (cf. |Mestel fc Spitzer (1956)| ; [Norman 



|fc Silk (19801 ; Falgarone fc Puget (19861 [McKee (19891) . Potential energy is liberated 



by overall cloud contraction and/or star formation (with winds); for quasi-equilibrium the 
respective rates would depend on how fast waves decay. At present, estimates of decay rates 
rely on analytic calculations of linear damping by ion-neutral friction, nonlinear single-wave 
steepening (see |Zweibel fc Josafatsson (1983)| and references therein), or simply dimensional 
analysis in terms of the internal velocity dispersion and size (see [Black (1987)| and references 
therein) . 

Previous theoretical studies of MHD waves in molecular clouds have concentrated on 
the linear, adiabatic, WKB limit in analyzing wave dynamics ( Fatuzzo fc Adams (1993" 



[McKee fc Zweibel (1995) ; Zweibel fc McKee (1995) ). This approach is accurate when the 



amplitude of the wave is small, when the wavelength of the waves is much smaller than 
the size of the system, and when the wave damping is weak. Since the first two of these 
assumptions do not strictly apply in molecular clouds, and the third may not either, it is 
interesting to see if we can relax them somewhat. In this paper, we undertake to study the 
full nonlinear development of moderate-amplitude MHD disturbances in a self-gravitating 
system, using numerical simulations. We concentrate our attention on the most basic 
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questions of how well clouds can be supported against gravity by nonlinear disturbances, 
and how long it takes MHD turbulence to decay. By employing simulations, we can 
generalize ideas of support by simple Alfven waves to include arbitrary self-consistent 
disturbances in the magnetic field, fluid flow, and density structure. We can test how far 
the linear-theory predictions for wave support carry over to the nonlinear regime, and go 
beyond the purview of simple linear theory to investigate the growth of structure, cascade 
of energy between scales, and associated dissipation. 

Expedience demands some sacrifices of realism for this first study. Our most severe 
simplification is to restrict the motions to plane-parallel geometry, so all dynamical variables 
are functions of one independent space variable and time. Thus we allow for transverse 
motions but not transverse derivatives. Another idealization which is less serious for the 
large-scale motions we consider is the neglect of ambipolar diffusion. We shall also assume 
the gas is isothermal. 

The plan of this paper is as follows. In §2 we review basic observational results 
for molecular clouds and the implied timescales in the context of simple theoretical 
considerations of cloud stability. We then present the equations we shall solve, our 
numerical method, and various tests used to verify its performance (§3). In §4 we describe 
the series of simulations we have performed. Finally, we apply our results to astronomical 
systems, discuss directions for future research, and summarize our conclusions (§5). 



2. Summary of Observations and Theory 

2.1. Cloud Properties and Scalings 

At present, the overall dynamical characterization of molecular clouds (e.g. delayed 
collapse, quasistatic contracting or expanding equilibrium) remains uncertain. In particular, 
the uncertainties about formation mechanisms and difficulty of assigning ages hampers 
efforts to deduce evolutionary trends (e.g. [Elmcgrccn (1993)]) . Arguments about the 



dynamical state are therefore indirect. The presence of several-million- year-old stars within 
cloud complexes indicates they have survived long enough to exceed the characteristic 
gravitational collapse times of at least the (higher-density) clumps; yet, a general internal 
collapse has not occurred. While cloud complex ages may not exceed the collapse times 
at their (lower) mean density, they are likely stabilized by something stronger than gas 
pressure because wholesale collapse on the dynamical times would lead to an unacceptably 
large galactic star formation rate (see review of |5hu et al. (1987)1) . Indeed, the presence 



of substructure at all scales within molecular clouds reveals that something much more 
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interesting than a nearly-pressure-free collapse is taking place. 

Spectral observations of molecular transitions show hypersonic widths everywhere 
except in the densest regions. These broad, often non-Gaussian molecular lines ( Falgarorii 
fe Phillips (1990]| ; [Miesch fc Scalo (1995)|) , give evidence that the gas they trace is in a 
highly turbulent state, with the turbulent amplitudes dependent on spatial scale. For cloud 
complexes, and massive component clumps and cores within them, internal linewidths 
a v typically increase with size scale R as a v (R) oc R 1 ^ 2 (this type of scaling was first 
pointed out by [Larson (1981)| ; see also Pame et al (1986)1 , ISolomon et al (1987)| , Falgaroiiq 1 
|fc Perault (1987)1 , |Myers fc Goodman (1988b)| ). Together with the typical observed 
scaling of mean density with size as n oc R^ 1 ( [Larson (1981)[ ), this relation implies that 
internal velocities within these structures are roughly virial, consistent with a response 
to self-gravitational confinement. The internal kinetic energy densities in the lower-mass, 
non-self-gravitating clumps (e.g. Carr (1987) 



Loren (1989a]j |Loren (1989"F)l [Herbertz 



Ungerechts, fc Winnewisser ( 1991 )| , Falgarone, Puget, fc Perault (1992)| , [Williams, Blitz, fc" 



Stark (1995)|) are instead comparable to the mean internal energy density of the surrounding 



GMC, suggesting a balance of stresses at the interface between the clump and an interclump 
medium - "pressure confinement" ( |Bertoldi fc McKee (1992)]) . The internal velocity and 
density structure in clouds points to the importance of magnetic fields; however, the lack of 
detailed information on field strengths and topology makes it difficult to decide how close to 
equilibrium any cloud, component clump, or core really is (see |Goodman fc Heiles (1994~)| ). 
Below we outline cloud observed properties, define relevant timescales necessary for scaling 
our numerical simulations, and discuss the corresponding simple stability requirements 
assuming uniform conditions. 

The average linear dimension L is typically 40pc for GMCs, lOpc for dark cloud 
complexes, and up to a few pc for the component clumps within these complexes that 
contain most of the mass. The volume-averaged density nn 2 is typically 25 — 100cm -3 in 
cloud complexes, and 10 3 cm -3 in clumps. Kinetic temperatures range from 10 — 50K. This 
implies an isothermal sound speed c s = JkT/jj, = 0.19 — 0.41km s" 1 , where /i = 2Am p 



(e.g. [Blitz (1993)| ; |Cernicharo (1991")]) . A convenient measure of the importance of magnetic 
fields is the parameter /3, defined by 



c 

v 2 a 



B 2 /4irp 



0.021 



( T 



V10K7 V10 2 cm 



n H2 



(—) 

V 10/iG J 



(3 < 1 implies a magnetically-dominated regime. Generally (3 ranges between 0.1 and 0.001 
in molecular clouds, reaching unity only in cloud cores. Notice that under isothermal, 



J Our (3 differs by a factor of 2 from the usual plasma (3 P = (gas pressure) / (magnetic pressure). 
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field- freezing conditions, (3 decreases as a cloud contracts homologously (B oc n 2 / 3 ). It is 
therefore natural to expect the importance of magnetic forces to increase relative to gas 
pressure gradients as an efficiently-radiating cloud condenses out of the interstellar medium. 

The sound wave crossing time over a scale L is 

L ( L \ / T \~ 1/2 

a characteristic gravitational collapse time is 

<c=f^y /2 = 9.9f^y 1/2 Myr, (3) 



GpJ ' V10 2 cm- 3 

and the Alfven wave crossing time over a distance L along a uniform field B is 



( * = ^r^f =7 - 6 feH»J UssJ Myr - (4) 



2.2. Stability Requirements in a Uniform Medium 

With the definition of collapse time in equation (|3|), the ratio of structure scale L to 
the Jeans length Lj = c s (rr / 'Gp) 1 ^ 2 satisfies L/Lj = t s /t c , so structures with t s > t c are 
above the Jeans limit and cannot be supported by thermal pressure alone. The importance 
of a mean field B to cloud dynamics can be measured by comparing the Alfven crossing 
time tA = L/va to the characteristic collapse time t c , 

(5) 

assuming a uniform field strength and density. When perturbations of wavelength L are 
applied to a cold medium with the wavevector perpendicular to the field, an analysis 
analogous to the classical Jeans treatment ( |Chandrasekhar fc Fermi (1953)| ) shows that 




for ^ < 1 the perturbation is stable when t\/t c < 1 and unstable otherwise. Thus for a 
uniform field strength Bq, any cross-field column of H less than 4 x 10 21 cm~ 2 (i?/10/iG) will 
be cross-field stable. This borderline-stable column is consistent with the observed column 
density seen for many clouds ( |Larson (1981) ) provided the field strength B w 25/iG. The 



cross-field stability criterion £aAc < 1 is equivalent to the requirement that the mass-to-flux 
ratio in a sphere carved out of this medium satisfy M/$ < 1/ (3G 1 / 2 ). 

For density perturbations with wavevectors parallel to the mean field B , the Jeans 
gravitational stability analysis including just thermal gas pressure is incomplete if, as has 



- 8- 



been proposed, Alfven waves provide the primary cloud support along the mean field. 
In quasilinear theory (cf. pewar (1970)| , [McKee fc Zweibel (1995)]) , we can estimate 



the importance of waves in opposing gravitational collapse along B by performing a 
Jeans- type analysis with k || B and pressure supplied by Alfven waves according to 
^Pwavc = -Pwave,o(p/po) 1/2 (assuming t A < t c , t s ). Here, P wavc ,o = |<5B| 2 /87r is the pressure in 
the fluctuating field, and we assume Alfven wavelengths short compared to the wavelength 
of the density disturbance. Defining the Jeans number rij = L/Lj, the criterion for stability 
on a scale L works out to be 

where the Alfven wave surface energy density is E w = J dx^(p\5\\ 2 + |<5B| 2 /47r) = 
L(|5B| 2 )/4tt. 

An extension of the above quasilinear theory to the nonlinear regime (Ew 3> pLc 2 ) 
would argue that for arbitrary rij, clouds could be stabilized against collapse along B 
whenever there is sufficient wave energy (in the absence of wave decay) . Properly, we do not 
expect the analysis above to extend to the nonlinear regime (indeed, while our nonlinear 
simulations do show correlations of wave pressure and density, they do not obey any simple 
law like P wave oc p 7p ). However, any argument that balances wave energy in the slab against 
gravitational potential energy - for example, a virial analysis - will yield the same scaling 
as equation ([|), since binding energy is proportional to njpLc 2 (see eq. [^]). The result of 



Pudritz (1990)| is the same as equation (^|) up to the factor 1/4. Since our simulations begin 
from uniform conditions, we find it convenient to use the quasilinear-theory prediction of 
equation (^) as a reference point. Of course, this "pseudo- Jeans" analysis, or any analysis 
in terms of an initial wave energy, becomes inapplicable if random motions dissipate rapidly 
compared to the collapse timescale. 



2.3. Wave Dissipation Mechanisms 



The dominant linear damping mechanism in molecular clouds is ambipolar diffusion 



see 



McKee et al. (1993)| and references therein). Ambipolar diffusion prevents propagation 
of Alfven waves with frequencies uja = k • va higher than a critical frequency uj c = 2u ni , 
where the neutral-ion collision frequency is v ni ~ rij 1.5 x 10 _9 cm 3 s _1 ( [Kulsrud fc Pearce 
1969) , Nakano (1984) ). Estimates of v n i indicate that waves can propagate at wavelengths 



well below O.Olpc flMyers fc Khersonsky ( 1995 )|) . Ambipolar diffusion also damps Alfven 
waves with oja < oj c at a rate E^/E^ = —u\/u ni . Writing the ionized fraction as 



^i/nH = Kn H , the ambipolar-diffusion damping time for wavelengths A is shorter than 
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the gravitational collapse time t c when 

x a a ( K Y 1/2 ( P Y 1/2 ( T \ 1/2 ( n Ha v 1 ' 2 

*< 14 Ln- B _-.- a /J 77777 77777 PC. (7) 



10- 5 cm- 3 / 2 / V 0.01 7 V10K/ V10 2 cm- 3 



For ionization principally by cosmic rays at the fiducial rate lO -17 * -1 , K ~ 10 -5 ( |McKee| 
|et al. (1993)|) , which for typical parameters can imply significant frictional damping on 
dynamical timescales for all wavelengths. Including the ionization produced by UV photons 
increases K and therefore decreases the frictional damping rate. For /3 < 1, fast MHD waves 
frictionally decay at a rate comparable to Alfven waves. Slow MHD (essentially acoustic) 
waves suffer little frictional decay, but may be subject to other linear dissipation such as 
radiative damping. 

Nonlinear (5B/B ~ 1) damping rates due to steepening of compressive (fast 
magnetosonic) and transverse (Alfvenic) isolated MHD wave trains have been compared by 
Zweibel fc Josafatsson (1983 )| . The former steepen at a rate kSv, while the latter steepen 



at a rate ^ k5v 2 /vA', transverse waves damp more slowly because pressure variations are 
second order in the wave amplitude rather than first order. Isolated wave trains deposit their 
energy in shocks in approximately a wave-steepening time. (Note that parallel-propagating 
slow MHD waves with 5v/c s ^> 1 do not exist; they shock immediately.) Because Alfven 
waves have the lowest nonlinear damping rates, it has been suggested that the observed 
supersonic linewidths owe their existence to these waves ( [Zweibel fc Josafatsson ( 1983 )|) . 

Isolated circularly polarized Alfven wave trains form a special case because they are 
exact solutions to the compressible ideal MHD equations and therefore suffer no nonlinear 
steepening. They are subject to a parametric instability, however, that causes them to 
excite compressive and magnetic disturbances over a range of uj and k QGoldstcin (1978] ). 



The growth rate of the parametric instability is of order the wave frequency when (3^1 
and the dimensionless amplitude of the wave is of order unity. In the small-amplitude, 
low-/3 limit this instability becomes the well-known decay instability ( [Sagdeev fc Galc~ 



1969)|) in which a forward-propagating circularly-polarized Alfven wave decays into a 
forward-propagating acoustic wave and a backward-propagating circularly-polarized Alfven 
wave. Circularly polarized Alfven waves cannot be relied upon, therefore, to provide slowly 
damped motions in molecular clouds. 

In truth, many waves are present in a molecular cloud simultaneously, and interactions 
among waves are as important as the steepening (essentially a self-interaction) of a single 
wave. Wave interactions transfer energy between scales; when energy is transferred to the 
smallest scales, dissipation occurs. The energy spectrum that develops depends of course 
on the allowed wave families in the fluid. The most familiar example of this process is the 
Kolmogorov cascade in an incompressible, unmagnetized fluid, which leads to an energy 
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spectrum dE{k)/dk oc k~ 5 ^ 3 . For magnetized, incompressible fluids (c s /t>A ^ 1), a theory 
of weak turbulence based on coupling of shear Alfven waves has recently been developed 
by [Sridhar fc Goldreich (1995)] (SG). The cascade to small scales occurs predominantly 
perpendicular to the mean magnetic field, but the weak theory becomes invalid at large k± 
due to increasing nonlinearity. A theory of critically-balanced strong turbulence, again for 
incompressible media, has been proposed by |Goldreich fc Sridhar (1995)1 (GS). 

Turbulence in molecular clouds lies in a rather different part of parameter space than 
that considered by Goldreich & Sridhar, so their theory cannot be directly applied. The 
SG+GS theory was developed for application to scintillation in the ionized interstellar 
medium, where c s /va ^ 1- In molecular clouds, by contrast, self-gravity and efficient 
cooling have conspired to produce a highly compressible medium (c s /v\ <C 1) in which 
turbulence is probably strong (Sv/va ~ 1) at the largest scales. In low-/? media, compressive 
disturbances are easily excited and damp rapidly via shocks. In general, a low-/? plasma can 
be expected to generate both an anisotropic shear- Alfven cascade, similar to that studied 
by SG+GS, and compressive disturbances and shocks (cf. Ghosh & Goldstein (1994)). 
Only with high-resolution two- or three- dimensional numerical studies of compressible 
MHD turbulence will it be possible to delineate the circumstances when one or the 
other mechanism provides the primary dissipation for Alfvenic disturbances in conditions 
appropriate for molecular clouds. In this first study, we have adopted a compressible 
equation of state but use the the simplest possible geometry (slab symmetry) that allows 
for transfer of Alfven -wave energy to compressive motions. Insofar as our model allows for 
compressibility but not a perpendicular Alfven -wave cascade, our study is complementary 
to that of Goldreich & Sridhar. 



3. Basic Equations, Numerical Method, and Tests 

In this paper we consider the simplest possible system in which nonlinear MHD waves 
and gravity can interact. We solve the equations of self-gravitating, compressible, ideal 
MHD: 

<9v , VP VB 2 (B-V)B „ ± 

— + ( v • V v = — + 1 ; - V</>, (8) 

at p Sirp Arrp 

^ = -V-(pv), (9) 

^ = Vx(vxB), (10) 
V 2 = 4vrGp. (11) 
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The equation of state is isothermal, 

P = c 2 sP , (12) 

with c s = 1 throughout. 

We assume "slab symmetry" , that is, all variables are a function of one independent 
spatial variable x and time t. Thus derivatives in the transverse (y and z) directions are 
zero, and the only derivatives that survive in equations (|8[) to (|TTJ) are in the longitudinal 
(x) direction. Transverse velocities v_|_ and magnetic fields Bj_ also survive and vary with x, 
making our scheme "1 + 2/3 D". Notice that because of the symmetry and the constraint 
V • B = we have B x (x,t) = const. The boundary conditions are periodic in a slab of 
thickness L — 1, i.e. a fluid element leaving the model at x — —L/2 reenters at x = L/2. 
For periodic boundary conditions, p — ► p — p on the right-hand-side of equation ([11]) (see 



Our numerical method is a one dimensional implementation of the ZEUS code described 
by |3tone fc Norman (1992ajl and Stone fc Norman ( 1992b) . The hydrodynamical portion 



of the code is a time-explicit, operator-split finite difference algorithm on a staggered mesh. 
Density and internal energy are zone-centered, while velocity components dwell on zone 
faces. The magnetic portion of the code uses the Method of Characteristics to evolve the 
transverse components of the magnetic fields and velocities in a manner that assures the 
successful propagation of Alfven waves. 

The gravitational acceleration, when used in our model, is calculated by taking 
the Fourier transform of the density and, for each component of the potential 4>k, 
setting (fik = 4:7iGpk(Ax) 2 / (2 cos(kAx) — 2), where Ax is the zone spacing. This 
gravitational kernel ensures that Poisson's equation is satisfied in finite-difference form. 
The gravitational acceleration is calculated by taking the inverse transform and setting 
9i = —((fii+i/2 — Ax. We have confirmed that selected spatial modes obey the Jeans 

dispersion relation as a test of the gravitational portion of the code. 

As a further comparison with linear theory, we have verified the scalings -P wave p 3 ^ 2 
and P W avc p 1 / 2 for the Alfven wave pressure in, respectively, an adiabatically contracting 
cloud and a wavetrain propagating along a density gradient QMcKee fc Zweibel (1995)|) . 



To verify the first scaling, we imposed a slow contraction of the coordinate system in a 
simulation containing a low-amplitude Alfven wave and found good agreement with the 
adiabatic wave amplitude-density relation. To verify the second scaling, we imposed an 
external gravitational potential to set up a density gradient and forced a low-amplitude 
Alfven wave at the bottom of the potential well. Again, we found that the wave amplitude 
for the outward-propagating wavetrain agreed well with the linear "polytropic" theory. 
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While our numerical method thus produces good agreement with linear theory, the 
experiments we present below involve the evolution of highly nonlinear MHD systems. It 
is hard to find good code tests for such systems, since there are few exact solutions known. 
Tests that exercise an MHD code's ability to handle the full family of MHD discontinuities 
have been described by [Ryu fc Jones (1995)1 ; we have performed these tests and obtained 
satisfactory results. 

Another direct and germane nonlinear test of the code is provided by the parametric 
instability of circularly polarized Alfven waves in a compressible fluid described by Goldstein 
|(1978)| . We have simulated the instability of circularly polarized Alfven waves and verified 
that the growing disturbances obey Goldstein's analytic dispersion relation. We have also 
simulated the steepening of an elliptically polarized Alfven wave and found good agreement 
with the calculations of |(Johen fc Kulsrud (1975)| . 

Finally, we have convergence tested the code for many of the simulations described 
in this paper. An example is shown in Figure 1, which displays the evolution of "wave" 
energy (kinetic and magnetic energy) in a suite of simulations that begin with a "random" 
set of transverse velocities and magnetic fields (see §4.1 for details of the initial conditions). 
The simulations vary only in their numerical resolution; the initial conditions are identical 
in the sense that their Fourier transforms are the same. Evidently as resolution increases 
the energy evolution converges. This is because numerical diffusion decreases as resolution 
increases, while the real sources of dissipation - isothermal shocks- become increasingly 
well resolved. Provided we use sufficient resolution (of order 512 zones in this case), the 
energy evolution will not depend qualitatively or quantitatively on numerical effects. 

4. Simulations 

We have performed several different types of simulations using the basic model 
described in the last section: a periodic, one- dimensional system with magnetic fields 
and an isothermal equation of state. This model may be thought of as mocking up a 
piece of a molecular cloud that is dominated by a mean field lying in the longitudinal 
(x) direction, although it also serves as perhaps the simplest possible context in which to 
study the interaction of Alfven waves and self-gravity. All the simulations begin with a 
uniform initial density p = 1, zero longitudinal velocity (v x = 0), and zero mean transverse 
field (By = B z = 0; we have verified that small mean transverse fields do not change 
our results qualitatively). Each simulation can be characterized by the (unchanging) 
strength of the mean longitudinal magnetic field B x , expressed in terms of the parameter 
f3 = cl/v\ = Cg/ (Bl/A-np) (see eq. |l|]; we include only the mean field component B x in 
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Fig. 1. — A convergence test. The figure shows the evolution of the sum of the magnetic 
and kinetic energy in simulations that begin with an initially random transverse field and 
velocity, but with different numerical resolution. This shows that the evolution does not 
depend on numerical effects for large enough resolution. 
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va)- The numerical resolution for the "standard" runs described below is 2048 zones; all 
other runs were done with 512 zones. Four different types of simulations are discussed in 
the subsections below. They are: the free decay of an initial wave spectrum (§3~T); a forced 



equilibrium in which the model is stirred at the same rate that it dissipates energy 
the free decay of an initial wave spectrum in the presence of self-gravity (§ |4.3|) ; and a forced 
equilibrium in the presence of self-gravity (§ |4.4j ). 



A useful diagnostic for the simulations described below is what we shall call the wave 
energy. It is defined as the sum of the kinetic energy and the energy of the perturbed 
magnetic field, E w = E K + E B , where 



r L/2 I 

/ dx-p(\v ± \ 2 + v 2 x ) (13) 

J-L/2 2 



r L / 2 IB i I 2 

E B = / dx 1 -^. (14) 



E K 

and 

Eb= I 

-L/2 87T 

Notice that E-&, E B , and E-w are in units of energy per unit surface area because of the 
slab geometry. For reference, the rms velocity and magnetic field perturbations satisfy 
(|5v/c s | 2 ) 1/2 = (2E K /pLc 2 )^ 2 and (\5B/ B x \ 2 ) 1/2 = ^' 2 {2E B / pLc 2 ) 1 ' 2 . The definition of p 
implies that v/v\ — f3 l ^ 2 v/c s for any v. Time is given in units of the sound-crossing-time 
t B = L/cs, and the Alfven -wave crossing time t& = L/v^ = t s f3 l l 2 . For self-gravitating 
simulations, the collapse time is t c = t s /nj, where rij = L/Lj (see eqs. [§], @, and ||). 



4.1. Free Decay 

First consider a non-self-gravitating model containing a spectrum of waves in the initial 
conditions that are allowed to freely decay thereafter. The "standard" run has (3 = 0.01 
and initial wave energy E w = 2E K = 2E B = lOOpLc 2 . With this initial energy, the initial 
field line distortions have a dimensionless amplitude (|(5B/^? a; | 2 ) 1//2 = 1. 

The initial transverse velocities and magnetic fields are a "random" superposition of 
parallel-propagating Alfven waves. They are drawn from a Gaussian random field with 
power spectrum (Iv^l 2 ) = (\'Bj_ i k\ 2 /47rp) oc k~ 2 for 2n/L < \k\ < 32(2ti/L). This power 
spectrum is special in that it approximates the natural power spectrum to which other 
initial power spectra decay, according to our numerical experiments, and also in that it 
corresponds to a velocity-dispersion/size relation a v (RY^ 2 oc R 1 ^ 2 in one dimension which 
is consistent with one of Larson's laws (cf. [Larson (1981)| ; §2.1| ). The power spectrum is 



steeply declining, so almost all the energy is at the largest scales. 
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Fig. 2. — Time evolution of the transverse magnetic + kinetic energy and the magnetic 
energy in the "standard" decay simulation. 
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The evolution of the total wave energy and the magnetic energy are shown in Figure 
2. The initial transverse magnetic fields give rise to a nonuniform magnetic pressure that 
is large compared to the gas pressure. This in turn produces longitudinal accelerations, in 
addition to the transverse accelerations caused by magnetic tension. The transverse kinetic 
and magnetic energies remain nearly in equipartition and oscillate out of phase, yielding a 
smoothly-evolving total wave energy; the longitudinal kinetic energy is only a few percent 
of the total. The wave energy declines to half its initial value by t = 0.2t s = 2.0t A . The 
transverse motions and fields lose their energy when the magnetic pressure due to the 
transverse fields does PdV work on the fluid, 

d f L / 2 J (I . |2 |BJ 2 \ fL/2 d f\B ± \ 2 \ 

dtL m dx \2 p ^ + ^rj = L /i dxv ^[-^r)^ (15) 

producing longitudinal motions that dissipate in shocks. Toward the end of the simulation 
almost all the wave energy is concentrated at the largest scales, thus giving the regular 
oscillations in transverse magnetic energy seen in Figure 2. 

Figure 3 shows density, longitudinal velocity v x , magnetic pressure, and the field line 
geometry of the simulation at t — 0.5t s = 5£a- There are large density contrasts present, 
with densities ranging over three orders of magnitude. Numerous shocks provide the 
dissipation. Notice that the magnetic pressure dominates the gas pressure (the gas pressure 
is numerically equal to the density, since cj. = 1). The field line structure in the x — i planes 
(i = y, z) is shown in the lower right panel of the figure, with the i field displacement Si 
defined by 

/X _g. 
dx — const., (16) 

where the constant is set so that Si = 0. Because most of the energy remains in the 
largest-scale Fourier components of the field, the field displacements are well ordered, with 
one maximum and one minimum. The field lines are nearly straight, and hence force-free, 
in between the kinks at density maxima. 

The variations in density at t — 0.5t s = 5£a can be further characterized by computing 
the fractional volume and fractional mass at p > p c (Figure 4). The figure shows that 10% 
of the mass (volume) is at p > 9.1p (2.7 p), 50% of the mass (volume) is at p > 3.1p (0.42p), 
and 90% of the mass (volume) is at p > 0.55p (0.012p). Thus the mass is concentrated in 
dense regions, while a significant fraction of the simulation volume is nearly empty. These 
results are typical of all our simulations, although density contrast increases significantly in 
the self-gravitating simulations discussed below. 

At time t = 0.5t s = 5t\, the best-fit slope s to the combined power spectrum 
|v±, fc | 2 + |B ±ife | 2 /47rp oc £r s yields s = 2.3 between k/(2n/L) = 1 and 100. The slope 
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Fig. 3. — Portrait of density p, longitudinal velocity v x , transverse magnetic pressure Bj_/8ir, 
and transverse field line displacements 5 y , 5 Z in the "standard" decay simulation at time 
t = 0.5t s = 5t A . 




Fig. 4.— Volume and mass fractions above density enhancement level p c /p, for the 
"standard" decay simulation at time t = 0.5t s = 5t\. 
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varies rapidly, however, with (s) = 2.4 ± 0.16, where () indicates an average over time. By 
smoothing the vj_ data over Gaussian windows of varying width and averaging over the box, 
and over time, we can compute a simulated "linewidth-size" relation (see §4.2). Averaged 
over time, the best fit relation is a v (R) oc R 7 . 

A vital question for the structure of molecular clouds is whether the wave energy 
concentrates in the regions of highest density. This can be answered by comparing |5B| 2 
with p after smoothing both on a scale A. In our standard decay simulation, magnetic 
pressure and density are weakly anticorrelated in a time-averaged sense, with slope 
| <5B | 2 oc p~ ai for essentially all smoothing scales A (in §4.3| we show that self-gravity changes 



this correlation). Thus magnetic pressure plays a role in confining the "clumps". But the 
clumps are not in equilibrium: they form, disperse, and accelerate, and material migrates 
from one clump to another. We do find that clump masses grow over time, however, as 
energy is transferred to larger scales and more coherent motions. 

We have surveyed free decay models with the same initial power spectrum over a 
variety of (3 and initial E\y. The evolution of the wave energy in a selected sample of these 
simulations is shown in Figure 5a. The initial dimensionless wave amplitudes range from 
(\5'B/B X \ 2 ) 1 ^ 2 = 1/4 to 8. Considering the subset of simulations where the initial rms field 
perturbations (j3Eyj/ pLc 2 ) 1 ^ 2 are 2 or less, for a given initial £7\y the decay is slower as the 
field strength increases (j3 decreases). For this same set of simulations, when we consider 
decays at a given (3, the fractional energy loss at t& increases as initial Eyj increases. 

Some of the simulations with initial wave energy of 400pLc 2 stop decaying near the 
end of their run. This is because they are able to make their way into a special dynamical 
state. In this special low-decay state, all the matter is concentrated into two narrow sheets 
that oscillate transversely. This special state is artificial, because its stability depends on 
the boundary conditions. We have tried duplicating the special state into a region of length 
2L, so that there are four narrow sheets, and introducing some small perturbations. We 
find that the system becomes unstable and decays to a new special state with only two 
sheets. The special state may also suffer higher dimensional instabilities which cannot be 
represented in the present simulations. While our simulations display a common trend 
toward concentration of energy and structure at the largest available scales, it seems 
unlikely that suitable conditions for the persistence of this state will be found in nature. 

It is not possible to define a "decay rate" from the free decay simulations because the 
decay rate depends both on the energy and on the internal state of the system. This is 
made clear in Figure 5b, which shows the instantaneous decay rate —E^/E^ for the same 
runs shown in Figure 5a. The heavy black squares indicate where the simulations begin; 
the decay rates quickly rise, then the energies and the decay rates fall. Except for initial 



-20- 




t / (L/c.) Wave Energy / pLe* 



Fig. 5. — (a) Time evolution of the magnetic+kinetic energy in a series of decay simulations 
with initial wave energy (transverse kinetic + transverse magnetic) of 25, 100, and 400, and 
initial f3 = c 2 Jv\ = 0.04,0.01, and 0.0025. The Alfven time is t A /t s = 0.05,0.1,0.2,0.4 for 
f3 = 0.0025,0.01,0.04,0.16, respectively, (b) The dissipation rate —E w /E w as a function 
of .Ew m the same simulations as (a). The simulations begin in the neighborhood of the 
solid black square. The dissipation rate initially increases, then decreases, as the simulations 
evolve and energy decays. Toward the end of many of these simulations the dissipation rate 
is sharply reduced because most of the power is in the lowest order {k = 2n/L) Fourier 
mode, which decays very slowly. The scale on the ordinate gives the decay rate in units i" 1 . 
For = 0.0025, 0.01, 0.04, 0.16, a decay rate t^ 1 = 20, 10, 5, 2.5t s " 1 , respectively. The straight 
lines show single wave steepening rates based on Cohen & Kulsrud (1974). 
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transients in a few of the high-/?, high-E'w cases, decay rates lie below t A . 

In Figure 5b, the straight lines show the nonlinear decay rate —Eyj/Eyj = 
27i(E w /pLc*)f3 1/2 t- 1 based on the single-wave steepening calculation of |Cohen fc Kulsrucfl 



1975)| , which has sometimes been used to estimate the dissipation rate within molecular 



clouds. This nonlinear steepening rate overestimates the decay rates we find by an order of 
magnitude or more, and does not display the same scaling as either the peak decay rates 
for different runs at a given j3, or the evolutionary decay rate for any given run. Indeed, one 
would not expect the Cohen & Kulsrud estimates to agree with our simulations, since: (a) 
the simulations contain a spectrum of waves, while the CK estimate is for a monochromatic 
wave train; (b) the disturbances in our simulations are highly nonlinear at the outset- they 
do not gradually steepen from a near-mode of the plasma. To the extent that a single decay 
rate 7 can characterize these freely-decaying systems as a function of E^/pLc 2 and /3, we 
note that the upper envelope of the decay rate in Figure 5b has better consistency with 
the scalings and approximate magnitude of the decay rate measured in the steadily forced 
simulations of the next subsection. 

The dissipation rate depends strongly on the initial power spectrum. We have 
performed a series of simulations with (v^l 2 = |B^ i / c | 2 /(47rp) oc k~ 2 for k m[n < k < k mauX 
and elsewhere, where k m i n /(2n/L) = 1,2,4,8, and 16 and k mSuX / (2tt / L) = 64. In all 
cases j3 = 0.01. The time evolution of the wave energy in these simulations is shown in 
Figure 6; evidently the smaller the scale at which the energy is injected the more rapid 
the dissipation. The peak decay rates are ~ 0.20;^. Our examination of the evolution 
of the density and field structure, and the power spectra of the transverse velocities and 
field (not shown), shows that while energy is lost overall, some fraction of the energy is 
transferred to modes with k smaller than the initial /c m j n . The final power spectra have 
l v _i_,A:| 2 , 1 Bj_ 5 fc 1 2 oc k~ 2 - 5±0 - 2 independent of where the cutoff k min / (2tt / L) was imposed in the 
input spectrum. We have also initiated decay simulations with alternative spectral slopes, 
and found that they too evolve toward spectra with slopes in the range —2 to —2.5. 



4.2. Forced Equilibrium Runs 

The MHD waves in molecular clouds may be trapped and preserved from the time of 
cloud formation, and the free decay simulations are designed to model that possibility. An 
alternative, however, is that they are continuously driven by processes inside the cloud, such 
as winds from young stars ( [Norman fc Silk (1980)|) . In this section we consider stochastic 
driving of transverse motions as a model for this process. 
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Fig. 6.— Time evolution of the wave energy (kinetic + transverse magnetic) in decay 
simulations where the initial energy power spectra extend between k m \ n / (2n/L) = 1, 2, 4, 8, 16 
(see labels) and k max /(27i/L) = 64. The figure demonstrates that decay proceeds more 
quickly when wave energy is stored in smaller-scale, higher-frequency oscillations. 
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The forcing is accomplished by introducing a pattern of transverse velocity perturbations 
Sv± at regular time intervals St. The velocity perturbations have the following properties: 
(1) |5v_|_ jfe | 2 oc fc 4 exp(— 4/c/fcpk), so that the power spectrum is peaked at k — k^ (2) the 
phases of the Fourier components of the perturbation are random; (3) the power spectrum 
is normalized and the mean component of the transverse velocities set so that 

f 1 

SE = J dx-p (Jv + Sv\ 2 - |v| 2 ) = const. (17) 



and 



J dxpSv = 0. 



The forcing interval St is typically set to be one one-hundredth of a sound crossing time, 
while k pk is 8(2ti/L) unless stated otherwise. 

The evolution of the standard forced run in shown in Figure 7, for which f3 = 0.01 and 
the forcing power SE/St = lOOpc^. The model equilibrates within a fraction of a sound 
crossing time. The steady state contains large density fluctuations, as in the free-decay 
runs. As in the free decay runs, the velocity and magnetic field power spectra evolve 
toward approximately oc k~ 2 (even though the model is forced at a scale well below the 
box size). The best fit slope for the combined power spectrum is s = 2.2 ± 0.12, averaged 
between times t/t s = 1 and 5. After saturation, the density power spectrum is flat below 
k = 8(2tt/L). 

We now turn to the "linewidth-size" relation in our simulations. In molecular clouds, 
the linewidth-size relation, one of "Larson's laws" ( [Larson ( 1981 )| ) , is a correlation with 



the approximate form a v (R) oc R 1 ^ 2 between the linear size R of a region and its internal 
velocity dispersion a v (R). Any linewidth-size relation is directly related to the power 
spectrum of the velocity: for \vk\ 2 oc k~ n ~ m in m dimensions, we have a v (R) oc R n ^ 2 
(assuming the emissivity of the gas is uniform). Larson's law is then a natural consequence 
of an n = 1 power spectrum. Power spectra with n = 1 occur in systems with an ensemble 
of velocity discontinuities, known as Burger's turbulence, and are observed in the solar 
wind ( |Burlaga fc Mish (1987)| ) and in simulations of supersonic hydrodynamic turbulence 



( [Passot, Pouquet, fc Woodward (1988)| , [Porter, Pouquet, fc Woodward ( 1992 )| ) . The power 



spectra in our simulations evolve to approximately |f_L,£;| 2 , |-B_i_,fc| 2 oc k~ 2 (independent of 
the forcing scale). By smoothing with a Gaussian window, we have directly verified the 
correspondence between the power spectrum and linewidth-size relation. For the standard 
run we find a linewidth-size relation of the form a v (R) oc R°- 57±om between scales L/512 
and L/2. 

We have performed a survey of forced runs to measure the saturation energy as a 
function of (3 and the forcing power. The results are shown in Figure 8; a fit to the results 
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Fig. 7. — Time evolution of the magnetic+kinetic energy and the magnetic energy in the 
"standard" forced run. The forcing is accomplished by adding random transverse velocities 
at fixed time intervals. 
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Fig. 8. — A contour plot of the saturation kinetic energy E K (denned as the average kinetic 
energy during the final two sound crossing times of the simulation) as a function of the 
forcing power and (3 = c 2 s /v\. The contours are located at E K /(pLc 2 ) = 2, 4, 8, 16, 32, 64. 
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that is better than a factor of 2 everywhere in the region surveyed is 



E * a67 



E w = 0.48 ^—3 1 P~"- Lb pLci. (19) 

We can use the results of Figure 8 to decide what input power is required to sustain a given 
level of turbulence. Solving for the required power based on the fit above gives 

E = 3 -°{M >m <2o) 

Since the model is in equilibrium the input power matches the dissipation rate. 

Notice that, apart from the factor of ft ' 24 ~ (c s /vaY^ 2 , the energy dissipation 
rate per unit volume is oc pa^/L, identical to what one would expect for homogeneous, 
incompressible turbulence. This is, in part, just dimensional analysis: a dissipation rate 
per unit volume must have units of a density times a speed cubed over a length. But it is 
not obvious which speed to use in completing this estimate, since there are three velocity 
scales in the problem: the sound speed, the Alfven speed, and the velocity dispersion. 
The simulations show that the correct speed is a power-law-mean of these three that is 
dominated by the velocity dispersion. 

An example of an alternative scaling for the Alfven -wave dissipation rate is the 
estimate of |Cohen fc Kulsrud (1975)] for the steepening of a single wave. The implied 



energy dissipation rate per unit area is E w ~ 2tt (E w / 'pLc 2 ) 2 j3 l / 2 pc^ s . The corresponding 
volumetric energy dissipation rate scales as (<j v /va) x (pa^/L), different from the results of 
our simulations. Thus the scaling of equation (|20|) is not a trivial prediction of quasilinear 
theory, even for compressible magnetized plasmas. 

The dissipation rate in the forced runs depends on the wavenumber at which the energy 
is injected (fc pk ). In the above runs the energy is injected near fc pk = 8(2n/L). A survey of 
the dissipation rate as a function of fc pk reveals that the saturation energy obeys 

E w oc fc p - k a3 °. (21) 

Together with equation(|T9|), this implies that 

E ex k° p f, (22) 

so that the dissipation rate increases as the forcing scale decreases. 

We can define a wave dissipation time tdiss i n equilibrium by taking the ratio of _E\y to 
E as given in equation (p0|). Including the forcing- wavelength dependence of equation (|22|), 
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we find 

u iss = ^ = o.87 ( ^V' 46 ( ^fV - 24 -. ( 

E w \L J \pLc 2 J c s 

We can compare the dissipation time for forced disturbances to the Alfven time t A 
or, for self-gravitating clouds, to the characteristic collapse time t c (eqs. ||, f|). Using 
E w /(pLc 2 s ) « 3a 2 /c 2 s , the former ratio is £ d i SS /*A ~ 0.5(A pk / 'L)^ 2 (v A / 'a v )(v A / 'c s ) 1/2 ; 
this ratio is generally large unless the forcing scale is very small. Similarly, 
^dissAc ~ 0.5(Apk/-^) 1 / 2 (wA/c s ) 1//2 (^jc s /cr t; ); note that a ratio smaller than unity does 
not necessitate collapse (see § |4.3|) . 



4.3. Self-Gravitating Decay Runs 

A crucial question for the evolution of molecular clouds is whether the nonlinear MHD 
waves interact with the self-gravitating compressive modes (here generically referred to as 
"Jeans modes" ) so as to prevent or delay collapse. We can examine this issue by turning on 
self-gravity in our simulations. 

Let us start with a few preliminary words about self-gravity in a periodic slab geometry. 
First, in one-dimensional, slab geometry an isothermal fluid always has a stable equilibrium, 
whereas in three dimensions a self-gravitating, isothermal fluid sphere is not generally 
stable (see jSpitzer (1968)| ). In non-periodic slab geometry, the equilibrium self-gravitating, 



isothermal density distribution is 

p = p sech 2 (z/z ), (24) 



where z a = c s j ^2nGp ( |Spitzer (1942)| ). In one dimension, the self-gravitating slab is the 
nonlinear outcome of the Jeans instability. Second, periodicity implies an important change 
in the physics of self-gravity: the adoption of the Jeans "swindle." When we solve Poisson's 
equation we make no allowance for the self-gravity of the background of material at the 
mean density. Thus we are really solving the equation 

V 2 = 4ttG(p - p) (25) 

rather than Poisson's equation. 

Consider a periodic slab geometry with initially uniform density. The self-gravity is 
such that there are nj Jeans lengths inside length L. If we assign a binding energy Eq = 
to the state with uniform density, it is possible to estimate the binding energy of the 
self-gravitating slab under the assumption that the density distribution is very nearly the 
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same as for a non-periodic self-gravitating slab, equation (|24]), and that there is only one 
slab present inside L. The result is 

E G , max ^1 - — - j pLcl (26) 

This is the equilibrium with the maximum possible binding energy for an isothermal fluid 
in one dimension. 

The energy evolution of the standard self-gravitating/decay run, with (3 = 0.01, initial 
wave energy E w = lOOpLcj., and Jeans number rij = 4 (so that E G)iriax ~ — 25pLc 2 ), 
is shown in Figure 9. | The background of transverse waves induces large density 
fluctuations within an Alfven crossing time, as in the non-self-gravitating case. These 
density fluctuations now have a gravitational binding energy, however. The density peaks 
grow and coalesce until at the end of the run (at t = t s = 10£a = 4t c ) there are only two 
peaks. As in the non-self-gravitating case, the field lines kink inside the density peaks and 
are nearly straight (i.e. force-free) in between. Since E G j E Gmax ~ 0.4 at the end of the 
simulation, the background of MHD waves evidently can delay collapse. 

A more detailed look at the evolution of the density structure in the standard 
simulation is provided in the right panel of Figure 10, which shows (p/p) 1 / 4 at intervals of 
£a = 0.1t s = 0.4i c . The density has been raised to the 1/4 power to reduce the contrast. 
The large density contrast induced by the the initial wave spectrum is evident at the first 
snapshot; the density maxima coalesce and grow over the course of the simulation. The 
left panel of Figure 10 shows a nonmagnetic simulation that begins with low amplitude 
white noise in p and v x . The nonmagnetic simulation provides a control for the magnetic 
simulation and shows that the presence of transverse MHD waves broaden the existing 
density maxima and delay or prevent their coalescence over several collapse times. 

How are the density and magnetic pressure correlated? As in § |4.1| , we have smoothed 
the density and magnetic pressure on a scale A, taken the logarithm of each, and measured 
the slope and strength of the correlation over a range of A. When smoothed on the smallest 
scales, the density and magnetic pressure are weakly correlated (averaging over time), with 
SB 2 oc p ' 1 . As A is increased, however, the strength and slope of the correlation increase. 
When A = L/5 the slope and strength of the correlation have approximately doubled. Thus 
on large scales the Jeans modes and Alfven waves do tend to interact in such a way as to 
resist collapse. Because the correlation is weak and variable, however, it is inappropriate to 
interpret it as an effective polytropic equation of state for the wave pressure. 



3 Since c s ,p, and L are fixed, the number of Jeans length is set by manipulating G so that, numerically, 
G = im 2 j. 
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Fig. 9. — Time evolution of the magnetic+kinetic, magnetic, and gravitational potential 
energy in the "standard" run with decay and gravity. There are four Jeans lengths inside 
the simulation. Thus t c = 0.25t s = 2.5t A . 
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Fig. 10. — Time evolution of the density structure in the "standard" self-gravitating run 
with decaying wave energy (right panel), and in an unmagnetized control run (left panel). To 
reduce the contrast, we display (p/p) 1 ^ 4 for both cases. Time intervals between the snapshot 
portraits are At — £a = 0.1t s = 0At c . The two dense structures which gradually coalesce in 
the magnetic simulation oscillate in the transverse direction. 
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Fig. 11. — A contour plot of the gravitational potential energy E G relative to the maximum 
gravitational potential energy for a fully-collapsed sheet -Eg, max (eq. 26), averaged between 
0.9 and 1.5 t c (the collapse time t c = t s /rij = 10t\/nj for (3 = 0.01, where nj = L/Lj is the 
number of Jeans lengths in the simulation) The ratio Eq/Eq^ max is measured as a function 
of the initial magnetic+ kinetic energy and nj. A large value of Eq/ £G,max corresponds to 
strong density concentration. The dashed curve shows the locus nj = (-Ew.init/ pLc^.) 1 / 2 /2 
(eq. [27]) predicted by the pseudo- Jeans analysis (§2.2) as the lower boundary of the collapsed 
region. See text for a discussion of the low-nj, high-i? w part of the diagram. 
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When can unforced MHD waves prevent gravitational collapse? To investigate this 
question, we have performed simulations with different initial energies and number of Jeans 
lengths. We measure the degree of collapse by comparing the time-averaged (from 0.9t c to 
1.5 £ c ) gravitational binding energy to the gravitational binding energy for a fully collapsed 
sheet, given by equation (p6|) . Figure 11 shows contours of this ratio {Eg I max) as a 
function of the initial wave energy and the number of Jeans lengths in the simulation. 

The maximum in Eg/ Ec y max in the upper left corner of Figure 11 is where collapse 
has occurred. This is qualitatively consistent with the pseudo- Jeans analysis of § |2.1| , which 
suggests that collapse should occur when the gravitational energy overwhelms the wave 
energy. Perhaps coincidentally, the pseudo- Jeans analysis is nearly quantitatively correct as 
well. The dashed line in the Figure shows the locus 



predicted by equation ([]), when Ey? ^> pLc^, as the lower limit in Jeans number for collapse 
to occur. This locus is indeed close to the boundary of the "collapsed" region. Other 
simulations show that the structure of this diagram is only weakly dependent on (3 over the 
range we have surveyed. 

There is also a maximum in Eq/ 'Eg, max m the lower right corner of Figure 11. In this 
region of parameter space the large initial wave motions drive large density fluctuations. 
These density fluctuations have an associated binding energy, which is large compared to 
the maximum possible binding energy E^max because E^max is small when rij is small. 
Thus, while the binding energy is near maximal, the system is not bound because the wave 
energy is comparatively large. 

The scaling with energy of the stability criterion (equation ||27|| ) and the dissipation 
timescale (equation [p^] ) leads to the following peculiar situation: in the stable region at the 
lower right of Figure 11, the dissipation time is shorter than the collapse timescale, while 
in the unstable region in the upper left of the Figure the dissipation time is longer than the 
collapse timescale. This apparent paradox is resolved when one realizes that stability is 
determined by the energy content of the cloud and not its energy loss rate. 



A final set of simulations considers stochastically forced runs with self-gravity. These 
runs are forced in the same manner as the nonself-gravitating, forced runs of ^4.2j , with the 




(27) 



4.4. Self-Gravitating Forced Runs 
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peak forcing at = 8 (2tt/L). The simulations are run without self-gravity for 0.5 sound 
crossing times, so that the model has a chance to equilibrate; then gravity is turned on. 

The standard run, with (3 = 0.01, nj = 4 (so E^max — — 25pLCg), and forcing power 
5E/5t = lOOpCg, is shown in Figure 12. The initial 0.5 sound crossing times is identical 
to the nonself-gravitating runs; collapse occurs as soon as self-gravity is turned on. The 
longitudinal motions associated with collapse cause a peak in kinetic energy near this time. 

Can collapse be prevented in these forced runs? We have performed simulations for a 
variety of nj and forcing powers and the results are displayed in Figure 13. The contours 
in Figure 13 lie at constant E G / E G ^ max . The pseudo- Jeans analysis leading to equation 
(|7|), together with the saturation energy predicted by the fit to the nonself-gravitating, 
forced runs (equation |TI|; the self-gravitating runs have a saturation wave energy that is 
comparable with this result, but generally smaller by a factor of < 2), predicts that collapse 
should occur for 

nj > 0.35/T 008 . ( 28 ) 

The predicted locus of marginal stability for f3 = 0.01 is shown as a dashed line in Figure 



13. The weak (3 dependence in equation (28) has been confirmed by runs not shown in the 
figure. The pseudo- Jeans analysis is thus a fair collapse predictor in this case as well. The 
corresponding lower limit for the power input needed to prevent collapse is 

E = 2A(3°-V J pcl (29) 

which for (3 = 0.01 has a coefficient 7.9. In any event, it is clear that transverse magnetic 
pressure can also prevent collapse in the forced runs. 



5. Discussion 
5.1. Applications 

In this section we translate our results into an astronomical context by applying them 
to a "typical" molecular cloud. This is a speculative venture, since our simulations employ 
slab symmetry and make a number of other drastic approximations. The exercise seems 
worthwhile, however, since the simulations are self-consistent and fully nonlinear, unlike 
earlier treatments of the problem. Are our results in harmony with the known, observed 
properties of molecular clouds? 

Consider a representative molecular cloud of mean linear dimension L = 20pc, number 
density 100cm -3 , kinetic temperature T = 20K, and line-of-sight velocity dispersion 
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Fig. 12. — Time evolution of the magnetic+kinetic, magnetic, and gravitational potential 
energy in the "standard" run with forcing and gravity. The strength of gravity is such that 
there are four Jeans lengths inside the simulation scale, and the forcing power is lOOpc^. 
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Fig. 13. — A contour plot of the gravitational potential energy Eq relative to the maximum 
gravitational potential energy for a fully-collapsed sheet -Eg, max (eq. [26]) measured at the 
end of the forced, self-gravitating simulations. E G / E G raax is found as a function of the 
forcing power and Jeans number nj = L/Lj. The dashed curve shows the locus estimated 
as the lower limit for collapse in equation (28). 
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<j v = 2km s _1 . For mean field strength B = 20/uG, (3 = 0.01. The sound-crossing, 

gravitational collapse, and Alfven -crossing timescales for this reference cloud are 

t s = 74Myr, t c = 9.9Myr, and t\ = 7.6Myr (assuming solar metallicity). The Jeans number 

= L/Lj = t s /t c = 7.5, and the ratio £aAc = 0.76. Thus, for a static uniform field, the 
cloud would be unstable parallel and stable perpendicular to the mean field (see § |2.2|) . If 
the projected area of the cloud is L 2 = 400pc 2 , then the kinetic energy per unit surface area 
is (3/2)pcr^L = 87pLc 2 . Assuming equipartition between kinetic and perturbed magnetic 
energy, the wave energy per unit surface area is twice the kinetic energy, Ew = !73pLc 2 . 
Finally, if the cloud's volume is L 3 , its total mass is 5.6 x 10 4 M Q . 

In the absence of self-gravity or any energy inputs, the cloud would, by assumption, 



evolve like the free decay simulations of §fO. In particular, it should evolve similarly to 



the "standard" decay of Fig. 2, which has Ew = lOOpLc 2 and (3 = 0.01, and is shown for 
comparison to the other decays as the central dotted curve in Fig. 5. To provide a more 
precise comparison, we have run a decay simulation starting from Ew = 1 73 pLc 2 s . After 
9.6Myr(= 0.16t s ), the wave energy would drop by factor « 2, so the velocity dispersion 
would drop by a factor of ^/2; after 41Myr by a factor ps 4, corresponding to a factor of 2 
drop in the velocity dispersion. At 15Myr, the cloud contains clumps. Fully 28% of the 
mass, but only 2% of the volume, lies at densities greater than 10 3 cm~ 3 = lOp, while 74% 
of the volume lies at densities below the mean. 

Now suppose the cloud's internal motions are forced (by winds from young stars, for 
example) at a scale A = L/8 = 2.5(L/20pc)pc, such that the mechanical power input V 
equals the wave dissipation rate C and the observed internal wave energy represents a 
saturated equilibrium. Again, we temporarily ignore self-gravity. We can use the forced, 
non-self-gravitating simulations to estimate the required input power and dissipation 
rate associated with a given equilibrium level of internal kinetic energy. In general the 
dissipation rate C should be roughly equal to the luminosity of the cloud in the important 
gas phase cooling lines, while the mechanical power input V would equal the sum of energy 
inputs from young stellar winds ( [Norman fc Silk (1980) ), Alfven waves emerging from 



collapsing, rotating cloud cores QGillis, Mestel, fc Paris (1979)| , |Gillis, Mestel, fc Paris 



|(1979)| , |Mouschovias fc Paleologou ( 1980 )|) , etc. In equilibrium V = C Using equations 



and (22), we find 



P - £ - 1 ^2knTi^) (lOto^)^) (m) (2^) L0 ' (30) 
A dissipation timescale can be estimated via tdi SS = E/E (eq. 

,3kms 




^ = 6.0(^-r) 1(^1 (ttt-) Myr. (31) 



-37- 



This dissipation timescale is comparable to the collapse timescale. 

In fact, our reference cloud is self-gravitating; it contains n 3 = (7.5) 3 = 410 Jeans 
masses. While computational expense prevents us from performing simulations with rij ^ 5, 
we find that equation ( p7|) is a good predictor of collapse on a dynamical timescale when 
nj < 5. Assuming equation fl27|) applies for more strongly self-gravitating clouds as well, 
we find that immediate contraction can be avoided when 

(|5B| 2 ) 1/2 > 18 ( -^-) ( -^r) /iG, (32) 
Xl 1 1 \ 20pci VIOOcm- 3 / P ' V ; 



corresponding to 



a v > 2.3 



\ 1/2 

J kms" 1 . (33) 



20pc J V 100cm" 3 , 

Our reference cloud with a v = 2km s _1 has velocity slightly below the threshold level, and 
in the absence of energy inputs would contract on a timescale ~ t c . 

The reference cloud could be supported in an indefinite equilibrium if it were supplied 
with energy at a scale A, as in the forced, self-gravitating simulations. Using equations ( p9[ ) 
and (|22|), we find that an input power of 




^ = 27 Ut£=5 77777 kr- L (34) 



or greater is required to support the cloud in a true equilibrium. A smaller supply of power 
would lead to gradual contraction. Notice that the power required to support the current 
level of turbulent energy in our reference cloud (eq. [3(]) would be insufficient to support it 
indefinitely against gravity. For the latter, the higher power level of equation ([34]) would be 
required, and the one-dimensional velocity dispersion would be raised to the level indicated 
in equation (|33|) . 



The "reference" cloud described above is meant to represent a cloud like Orion A 
( |Bally et al (1987)| ) or the Rosette nebula ( |Williams, Blitz, fc Stark (1995)| ). For clouds 
as large as this, we have had to extrapolate our fits to larger values of the Jeans number 
and input power than we were able to compute directly. On the other hand, smaller 
self-gravitating clouds such as Taurus-Auriga ( |Uernicharo (199IJD or Ophiuchus flLoren 



(1989a)|, [Loren (1989b)[) are within the parameter range we have simulated, and therefore 



the results presented in the figures of §4 can be used directly. Diffuse, high-latitude clouds 
may be best represented by our unforced, non- self-gravitating simulations. 

The power requirements for sustaining turbulence and counteracting gravity of 
equations (BO) and (B3) are reasonable for GMCs. For example, the total hydrodynamic 
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outflow momentum observed in Orion A (which has properties comparable to our reference 
cloud) is estimated at 32OM kms _1 ( [Fukui et al (1993~)| ), which for mean outflow velocity 



30km s 1 and outflow lobe size 0.75pc (cf. [b'ukui et al (1989~)| ) implies a characteristic power 
input of 33 L Q . Thus our results are in accord with observations. 

Finally, let us suppose that clouds evolve in quasi-equilibrium, and support can 
ultimately be ascribed to the power originating in young stellar outflows. Then by equating 
the required power V to the total wind mechanical luminosity GM.*M*/ R c we can obtain 
a total star formation rate A4 in the cloud. Here we absorb the details of the wind 
acceleration mechanism into the characteristic radius R c where the wind originates (see, 
e.g., |Shu et al (1994)|), and the uncertainties associated with the coupling of the wind to the 



rest of the cloud into the scale A. Scaling R c to 0.1 AU and using equation (PW), we find 



a. = 2xio- (j-x r^r (A) 1 " y-r r*r 



20pc/ VIOOcm- 3 / 10.01/ V 2.5pc / V0.1AU/V M © 



(35) 

This implies a star formation timescale M(cloud)/Ai* — 3 x 10 9 yr, which is about 300 
times the dynamical collapse time. It also implies that over a cloud lifetime of perhaps 
3 x 10 7 yr the cloud will turn about 1% of its mass into stars. 



5.2. Perspective 

It behooves us to remind the reader of some of the shortcomings of our treatment. The 
most severe is the use of slab (one dimensional) symmetry. In higher dimensions, dissipation 
rates may rise because new, transverse decay modes will be available. On the other hand, 
dissipation rates could fall somewhat because clumps, once formed, need not collide with 
each other, as they do in one dimension. Fortunately, an extension to two dimensions is 
immediately practicable. 

Although the use of slab symmetry is likely the leading source of error in our 
calculation, we have made other approximations that contribute as well, and that could 
be relaxed in future treatments of this problem. In this work, we have used an isothermal 
equation of state. Numerical work has shown, however, that realistic molecular cooling can 
have a significant effect on fragmentation during gravitational collapse (e.g. [Monaghan fc| 
Lattanzio (1991)|) . We also neglect ambipolar diffusion. This is likely to have particularly 



significant effects at small scales in clouds, where MHD waves cannot propagate. Even for 
Alfven waves on the scale of the GMCs, the ambipolar-diffusion damping timescale can be 
smaller than the dynamical timescale if the only source of ionization is cosmic rays at a rate 
10~ 17 s _1 . Overall, inclusion of ambipolar diffusion would tend to increase total dissipation 
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rates, and to steepen the power spectrum. 

Another idealization in our simulations is the periodic boundary conditions. These 
boundaries prevent any wave energy losses from our model clouds by radiation to an 
external medium, as would reflecting boundary conditions. Estimates indicate that under 
certain conditions wave radiation may dominate other energy losses from molecular clouds 
(cf. [Elmegreen [T985J ). For a linear amplitude wavetrain, the transmission coefficient at 
the cloud edge is small when the density changes sharply over a distance smaller than the 
wavelength. Since most of the energy in our simulations is at the largest scales, our results 
could be sensibly applied to clouds with edge gradients sharper than the inverse of the 
cloud size. 

Finally, we have neglected cosmic ray pressure and transport, radiative transfer effects, 
and, of course, feedback from star formation. It is not likely that these effects can be 
incorporated in any realistic way in the near future. Our forcing algorithm does model 
power input by, e.g., stellar winds, but only in a crude fashion. 

Nevertheless, our treatment does represent significant progress. As far as we are 
aware, it is the first fully self-consistent treatment of a turbulent, magnetically dominated, 
compressible, self-gravitating fluid. Our self-consistent simulations with realistic field 
strengths and velocities lead to a highly inhomogeneous state characterized by MHD 
discontinuities. This state bears little resemblance to any regime that has been well studied 
in the past. In particular, while some results carry over from quasi-linear theories and 
theories of incompressible, eddy-dominated turbulence, neither can adequately represent 
the internal dynamics of molecular clouds. 

5.3. Summary 

In this paper, we set out to explore the hypothesis that magnetic forces are crucial 
to the internal dynamics of Galactic molecular clouds. We were motivated by the 
proposal that the fluctuating velocity field in MHD waves is responsible for the observed 
hypersonic turbulence in molecular clouds, while the associated magnetic field fluctuations 
provide a pressure vital in supporting clouds against gravitational collapse and confining 
non-self-gravitating clumps. In particular, we were interested in whether observed molecular 
clouds represent dynamical equilibria, and more generally what is required to sustain a given 
level of MHD turbulence in the face of nonlinear dissipation and self-gravity. A quasistatic 
equilibrium is possible even with wave decay; cloud turbulence may be replenished for 
example by the gravitational potential energy liberated from rotating cloud cores and disks 
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when they collapse and accrete to form stars. 

This work uses numerical simulations to investigate large- amplitude (Sv/va ~ 1) 
MHD turbulence under density, temperature, and magnetic field conditions appropriate 
for Galactic molecular clouds. Computational expense has imposed some sacrifices of 
realism in this first self-consistent treatment of turbulence in a highly compressible, 
magnetically-dominated (c s /t>A <C 1) fluid. Our main idealizations are to restrict the 
dynamics to plane-parallel geometry with periodic boundary conditions, and to ignore 
ion-neutral slippage. This is the simplest possible model that incorporates magnetic fields 
and gravity self-consistently. 



We have performed four types of numerical simulations. In the first (§ [4.1[) , we evolved 
a spectrum of large-amplitude Alfven -wave turbulence to determine the free decay rate. 
In these simulations, we also investigated the spectral evolution of the turbulence and 
formation of density structure. In our second set of simulations (§ |4.2|) , we set out to 
evaluate the level of internal stochastic forcing (intended to model power inputs like young 
stellar winds) required to sustain a given level of turbulent kinetic energy. In our third 
set of simulations (§4.3|), we included self-gravity in model clouds initiated with differing 
levels of wave energy to establish collapse thresholds while allowing for turbulent decay. In 
our final set of simulations (j j4.4|) , we applied stochastic internal forcing to self-gravitating 



clouds to evaluate the power input needed to prevent collapse. As a template for using the 
dimensionless results of §|J we translate into physical units for a "reference" cloud in § [5T 



Our major conclusions follows: (1) We have confirmed that nonlinear 

disturbances in ideal MHD can support a model cloud against gravitational collapse, 
provided that the perturbed magnetic energy is maintained at a level exceeding the binding 
energy of the cloud (cf. eq.|3"2]]). Gravity is opposed by a gradient in the time- varying 



magnetic pressure due to the components of the field perpendicular to the mean field. 
(2) We have characterized the dynamical state of highly nonlinear, c 2 s /v\ <C 1 MHD 
systems. Such systems contain strong density contrasts, with much of the volume effectively 
evacuated and most of the mass concentrated into small regions. High density "clumps" 
form and disperse over time, with a secular trend toward increasingly large concentrations 
and coherent motions as energy cascades to larger scales. These nonlinear, magnetically 
dominated systems contain numerous MHD discontinuities, which naturally give rise to 
a wave energy power spectrum approximately |-Bj_,fc| 2 oc |f_i_,fc| 2 oc k~ s with s ~ 2, or 
linewidth-size relation approximately a v (R) oc R 1 ^ 2 . The shape of the power spectra at 
late times are essentially independent of the initial input spectrum or energy injection 
scale. (3) We have calculated a decay rate for nonlinear MHD waves in slab symmetry by 
equating forcing and dissipation rates in a saturated equilibrium. The dissipation time, 
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given in physical units in equation (pTj) , is longer than some naive estimates. In particular, 
when the cloud is stirred at scales comparable to its size, the dissipation time exceeds 
the "eddy turnover time" L/a v by {va/c s ) 1 I 2 times an order-unity factor, and the Alfven 
-wave crossing time by an additional factor v^/a v (see eq.|2"3fl). Peak dissipation rates for 
free decay obey approximately the same scaling in a v , va, and c s as the dissipation rate in 
saturated equilibrium. Because of the present restricted geometry and negligible friction, 
our computed dissipation rates are likely lower limits. 

We are especially grateful to Jim Stone for helping to initiate this project. We 
would also like to thank Bruce Elmegreen, Chris McKee, Phil Myers, and Phil Solomon, 
for interesting discussions, Alyssa Goodman, Jeremy Goodman, and Jerry Ostriker for 
comments on a draft of this paper, and the referee Ellen Zweibel for a discerning review. 
This work was supported in part by NASA grant NAG 52837. 
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